# Codigo para calculo da deflexao de membrana elastica circular
from dolfin import *

# Set pressure function
T = 10.0 # tension
A = 1.0  # pressure amplitude
R = 0.3  # radius of domain
theta = 0.2
x0 = 0.6*R*cos(theta)
y0 = 0.6*R*sin(theta)
#sigma = 0.025
sigma = 50 # verification
print 'Type the value for sigma(the width of the pressure function): '
#cin >> sigma

pressure = '4*exp(-0.5*(pow((%g*x[0] - %g)/%g, 2)) '\
           '      -0.5*(pow((%g*x[1] - %g)/%g, 2)))' % \
           (R, x0, sigma, R, y0, sigma)

n = 100 # approx. no of elements in radial direction
mesh = UnitCircle(n)
V = FunctionSpace(mesh, 'CG', 1)

# Define boundary condition w=0

def boundary(x, on_boundary):
    return on_boundary

bc = DirichletBC(V, Constant(0.0), boundary)

# Define variational problem
v = TestFunction(V)
w = TrialFunction(V)
p = Expression(pressure)
a = inner(grad(w), grad(v))*dx
L = v*p*dx

# Compute solution
problem = VariationalProblem(a, L, bc)
w = problem.solve()

# Plot solution and mesh
plot(mesh, title='Mesh over scaled domain')
plot(w, title='Scaled deflection')
p = interpolate(p,V)
plot(p, title='Scaled pressure')

# Find maximum real deflection
max_w = w.vector().array().max()
max_D = A*max_w/(8*pi*sigma*T)
print 'Maximum real deflection is', max_D

# Verification for "flat" pressure (big sigma)
if sigma >= 50:
    w_exact = Expression('1 - x[0]*x[0] - x[1]*x[1]')
    w_e = interpolate(w_exact, V)
    w_e_array = w_e.vector().array()
    w_array = w.vector().array()
    diff_array = abs(w_e_array - w_array)
    print 'Verification of the solution, max difference is %.4E' % \
          diff_array.max()

    # Create finite element field over V and fill with error values
    difference = Function(V)
    difference.vector()[:] = diff_array
    plot(difference, title='Error field for sigma=%g' % sigma)

# Calculate elastic energy in the membrane
energy = 0.5*inner(T*grad(w), T*grad(w))*dx
E = assemble(energy, mesh=mesh)
print 'Elastic energy of the membrane: %g' % E

# Hold plot
interactive()


